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We discuss the creeping motion of plugs of negligible viscosity in rough capillary tubes 
filled with carrier fluids. This extends Bretherton's research work on the in finite-length 



( Bretherton 1 961) 



bubble motion in a cylindrical or smooth tube for small capillary numbers Ca 
We first derive the asymptotic dependence of the plug speed on the finite length in the 
smooth tube case. This dependence on length is exponentially small, with a decay length 
much shorter than the tube radius R. Then we discuss the effect of azimuthal roughness 
of the tube on the plug speed. The tube roughness leads to an unbalanced capillary pres- 
sure and a carrier fluid flux in the azimuthal plane. This flux controls the relaxation of 
the plug shape to its infinite-length limit. For long-wavelength roughness, we find that 
the above decay length is much longer in the rough tube, and even becomes comparable 
to the tube radius R in some cases. This implies a much-enhanced dependence of the 
plug speed on t he plug length. Th is mechanism may explain the catch-up effect seen 
experimentally | Jlsmagilov k, Ying ) . 
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1. Introduction 



The field of multiphase microfiuidics has expanded greatly in the l ast decade because of 



(jGunther fe Jensen 20061 ) 



the numerous applications of the technique in chemistry and biology 
Transporting reactants through microchannels in the form of bubbles or droplets has 
many advantages over the traditional single-phase syste ms. These include enhance d mix- 



ing rate, reduced dispersion and higher interface areas | (jGunther fc Jensen 20061) . How- 
ever there are still many technological problems that remain unsolved. A particular prob- 
lem that aroused our attention is the control of spacing within a train of long droplets 
or plugs carried along the channel by a wetting fluid. Such plugs are separated from the 
tube by a thin lu bricating film, whose thick ness is controlled b y the speed of motion 



(jBretherton 19611) . It is commonly observed | (jlsmagilov fe Yina ) that the separation be- 
tween the plugs changes as they move. Eventually one plug can coalesce with its neighbor, 
a process known as the catch-up effect. To our knowledge there is no widely accepted 



explanation of this effect. 



A change in the distance between neighboring plugs is possible only if the fluid flux in 
the lubricating layer is different in the two plugs. Thus in order to understand the reason 
for the catch-up effect one needs to analyze the shape and dynamics of this lubricating 
layer. There are many physical mechanisms that could be responsible for the structure 
of the lubricating layer. In this work we analyze only one of them: irregularities of the 
channel shape. We argue that these irregularities, i.e. roughness of the tube wall, could 
be responsible for the fluctuations of the distance between the plugs. The influence of the 
substrate geometry or roughness on the lubricating layer thickness has been recognized 
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for several decades and has b een studied quite extensively (IStillwagon fc Larson 



Schwartz fc Weidner (1995) 



( 20011) ; 



Howell 



Kalliadasis. Bielarz fe Homsv 



2000); 



1988); 



Mazouchi fc Homsv 



2003)). In the present work we consider an effect beyond the direct affect 
on speeds. We find that the roughness enhances dramatically the sensitivity of the thin 
film width to the plug length, and thus produces speed variations between individual 
plugs that might eventually lead to catch-up coalescence events. 

The earliest discussion of the roughness effect on t he plug motion inside capillary 



tubes dates back to the seminal paper of Bretherton | (jBretherton 19611 ) . He suggested 
that the roughness can be important if its amplitude is comparable to or larger than 
the lubricating layer thickness. Therefore the roughness effect is usually negligible for 
large enough tubes whe re other physical mechan isms are more important (see e.g. 



(jGunther fc Jensen 20061 ) , | (jAiaev fc Homsv 20061 ) for review) . However as modern mi 



crofluidic channels become smaller, surface roughness becomes more important relative 
to the channel width. Probably the first experimental ob servation th at is attributed to 



the effect of roughness is reported in the work of Chen| ([Chen 19861) . where deviations 



from the Brcthcrton's predictions were observed for small capillary numbers Ca. 



On the theoretical side the most widely studied effect of the surface roughness corre- 
sponds to the limit where the roughness amplitude is significantly smaller than the lubri- 
cating film thickness. In this case the roughness effect can be accounted for via effective 
boundary conditions: the rough surface can be modeled as a smooth one wi t h non - 



zero slip length determined by the rou ghness properties ([Einzel. Panzer fc Liu 



Miksis fc Davisl(|1994l ); 



DeGennesI ([20021 ) ) . Krechetnikov and Homsy 



1990) 



( Krechetnikov fc Homsv 20051 ) 



have analyzed the effect of the slip length on the lubricating film profile, and showed that 
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the effect is negligible for the slip length smaller than the Bretherton thin film thickness, 
but leads to a dramatic increase of the lubricating film thickness in the opposite limit: 
in this case the lubricating film thickness becomes comparable to the slip length and is 
independent of Ca. In a broader context there have been a numerous studies on how 



the ch a nnel ge o metry affects the dynamics and shape of plugs. Wong et al. ([Wong 1 



(jl995h : 



Wong' 



(|1995| )) considered the plug motion in polygonal channels and showed 
that the shape of the plug strongly depends on the capillary number and that in the 
limit Ca «C 1, the fluid flux in the lubricating film does not depend on Ca. Hazel and 



Heil (|Hazel fc Heil 20021 ) confirmed this effect by numerical simulations, and showed that 
a similar behavior is observed in tubes with elliptical cross-sections as well. It should also 
be noted that the surface roughness can have a non -hydrodynamic effect on the p l ug mo - 



tion by changing wetting^ prop erties of the surface (jWenzell (|1936f ) 



Herminghausl (|2000h 



Bico. Tordeux fc Querd (|200ll) ) or dynamical destabilization of the lubricating layer due 



to fluctuations. 



In this work we consider capillary tubes with roughness smaller than the lubricating 
film thickness. In this regime the roughness cannot produce any significant effect on the 
width of the lubricating layer. However, as we will show, in the rough tube the decay 
length of the plug shape to its infinite-length limit is increased dramatically. In some 
cases, the decay length becomes comparable to or even larger than the tube radius. The 
ripples of the plug shape propagate from the front region to the rear region of the plug, 
leading to a much-enhanced dependence of the fluid flux on the plug length. On the level 
of the plug train this effect translates the variations of plug speeds to the fluctuations of 
the distance between them. Thus it can be responsible for the catch-up events. We also 
note that this is essentially a non-local effect. To our knowledge this effect has not been 
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discussed perviously. 
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We will illustrate this effect by using a very simple system, a single air plug of finite 
length moving through an almost cylindrical tube. For the sake of simplicity, through- 
out the paper we neglect the effects of gravity, inertia and surfactant concentrations, 
assuming that the only relevant stresses are due to surface tension and viscosity. We 
also limit our analysis to the small capillary number regime: Ca -C 1. We assume that 
the tube roughness can be modeled as weakly non-circular cross sections of the tube, 
which does not change in the longitudinal direction and can be described by a roughness 
function e(9) with 9 being the azimuthal angle. We will consider the situation where 
the deformation is smooth enough for the lubrication approximation to be valid, that is, 
e(6) <C R ■ Ca 2 / 3 , i.e. the deformation is small compared to the Bretherton's thin film 
thickness. However, we expect that our scaling results can be safely extrapolated to the 
regime: e(0) ~ R ■ Ca 2 / 3 . 



The structure of this paper is as follows. We start our analysis by revisiting the classi- 
cal Brcthcrton solution for the semi-infinite air plug mo ying in the per fectly cylindrical 



(Tcle tzke 1983 ) and extend it to 



tube. Then we follow the finite-plug analysis of Teletzke 
a wider range of plug length. The main finding of this section is the exponential suppres- 
sion of the speed corrections due to the plug length variations: 5U/U oc exp [—L p / (2 Loo)], 
where L n measures the length of the plug, and Loo = 0M3R ■ {iCa) 1 / 3 is the charac- 
teristic decay length of the plug shape to its infinite-length limit for the smooth tube 
case. In Section 3 we explain that the presence of small tube roughness produces a whole 
spectrum of relaxation modes with different decay lengths Aj. The relaxation mode with 
the decay length \ max ~ L p provides the largest contribution to the dependence on the 
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plug length. In Section 4 we estimate of the amplitude of the plug shape distortion due 
to the roughness. Then we estimate the effect of tube roughness on the finite plug speed 
in Section 5. We conclude in Section 6 by discussing the limitations of the model and 
proposing future research directions. 



2. Plug speed dependence on the plug length in smooth tubes 

Suppose we have a smooth cylindrical capillary tube with radius R. It is initially filled 
with a liquid of viscosity fx, which is called the carrier fluid in the following discussion. 
Now an air plug with negligible viscosity is forced into the tube and moves slowly together 
with the carrier fluid due to a pressure gradient along the tube. As discussed in detail 
later, the plug speed U is different from the speed V of the carrier fluid far away from the 



plug | ([Hodges. Jensen fc Rallison 2 004). In this section, we will investigate the effect of 
the plug length on the plug speed U. Both gravitational and inertial effects are assumed 
to be negligible. 



2.1. Bretherton-Teletzke mechanism 



In Bretherton's seminal work| (jBretherton 19611) . he studied the limit case of semi- 



infinite air plug motion. He showed that a dimensionless number Ca, called capillary 
number, controls the plug motion: 



Ca 



7 



(2.1) 



where 7 is an interfacial tension between the plug and the carrier fluid. Bretherton's 
steady state solution of a semi-infinite plug is illustrated in FigJTJ For convenience, we 
have chosen the plug as the frame of reference. In this frame, the tube wall moves with 



Plugs in rough capillary tubes 



7 




Figure 1. Sketch of a plug identifying Bretherton's quantities used in the text. 

velocity —U. The shape of a semi-infinite plug consists of three regions: a spherical cap 
(static region) is connected to a thin film region of almost constant liquid film thickness 
hoo via a transition remobilization region. Since the internal viscosity of the plug is zero, 
there can be no shear stress at the plug surface. Thus the carrier fluid in the thin film 
region moves without shear at the wall speed U. Bretherton also found that the lubri- 
catin g film thickness d ecays exponentially with a decay length hoc (3Ca) -1 / 3 , denoted as 



(jBretherton 19611) . Thus hoc is a measure of the length of the remobilization region. 
Both capillary force and viscous force are important in the remobilization region, where 
the plug shape is determined by an equilibrium between these two forces. In the limit of 
small capill ary numbers, Ca < 0.005 <C 1, Bretherton derived the dependence of hoc on 



Ca and R\ (jBretherton 19611 ): 

hoc = 0.643 R- (3Ca) 2/3 . (2.2) 
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Obviously, both L M and hoo are much less than R for small capillary numbers. Relative 
motion between the plug and the carrier fluid is governed by an incompressibility con- 
straint. This constraint states that, in the steady state, the total volume of the plug and 
the carrier fluid between two arbitrary cross sections along the tube doesn't change over 
time. Mathematically, it can be translated into the equation: ttR 2 V — ir(R — ft-oo) 2 U. 
By keeping only the lowest order approximation, we get: 

U/V ~ 1 + 2h oa /R 

= 1 + 1.286 (3C*a) 2/3 . (2.3) 

Strictly speaking, Eqn. (|2.3l) is only valid in the limit of the sem i-infinite plug. For a finite 
plug, we expect Eqn.fllS) to be modified as| (jTeletzke 19831 ): U/V = 1 + 2h*/R, with 
h* different from its limiting value hoo. Physically speaking, — h* represents the average 
carrier fluid flux per unit circumference. As shown in detail later, h* depends on the plug 
length and approaches hoo as we increase the length. 



The motion o f a finite-length air plug in a capillary tube was first studied by Teletzke 



( Teletzke 1 983). In the following discussion length, velocity, and pressure are made di- 
mcnsionless by using tube radius R, plug speed U, and capillary pressure j/R. Capillary 
number was kept small in Teletzke's analysis: Co < 1. By imposing a lubrication approx- 
imation and appropriate boundary conditions (non-slip condition at the tube wall and 
stress-free conditi on at the plug surface), the balance of viscous and capillary pressure 



gradients yields | ([Teletzke 19831 ): 

d 3 h 3Ca ■ (h - h*) 



(2.4) 



dz 3 h 3 

wherever \dh/dz\ <C 1. For small Ca, this includes both the thin film region and the 
remobilization region. For a plug of fixed length, h* is a self-consistent parameter whose 
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Figure 2. Quantities characterizing a plug of finite length, 
value is determined by the fact that remobilization regions must match with spherical 
caps at the plug ends. Suppose the front and rear remobilization regions are separated 
by a distance L p for a finite plug. Teletzke proposed a method to solv e the above dif - 



ferential equation numerically in the limit of short plugs: L p < lOL^ \ (jTeletzke 19831 ). 
Numerical results obtained by using Teletzke's method are plotted in Figj3] As expected, 
h* approaches ft.oo as the plug length increases. However, the trend is not monotonia 
The origin of this oscillating behavior will be clear in the next section. 



Teletzke's method did not address the asymptotic behavior for finite plugs with length: 
L p Loo. So the knowledge of the dependence of the plug speed on the plug length 
is incomplete at the current stage. We only understand two limit cases: Bretherton's 
semi- infinite plugs and Teletzke's short plugs with L p < lOLoo- In the next section, we 
complete this previously unstudied part of the plug speed problem, by determining the 



10 



1.3r 



1.2- 



1.1 



h*/h a 



Q. Zhang, K. S. Turitsyn, and T. A. Witten 



1.0 



0.9 



• ' 







Figure 3. Dependence of h* on the plug length L p obtained by using Teletzke's method, 
asymptotic behavior of h* for plugs with large L p . 



2.2. Asymptotic corrections to the semi-infinite plug speed 

Due to rotational symmetry, the problem of determining plug shape is reduced to 
calculating a lubricating film thickness . The governing equation of h(z) is given by 
Eqn. (|2.4|) . For small capillary numbers Ca, this differential equation is valid in both 
the thin film and the remobilization region, as noted above. Its solution has to match 
spherical caps at the ends of remobilization regions. By using the rescalings: h = h* ■ r\ 
and z — L ■ £, we may rewrite Eqn. (|2.4[ ) as: 

fl = H_ 1. (2.5) 

dt; 3 T] 3 

In the above substitutions, L is the remobilization length of the finite plug defined as: 
L = h* (3Ca) -1 / 3 . Comparing with Eqn. (|2.2p . we observe that r\ measures thickness in 
units comparable to the asymptotic thin film thickness h^,. Likewise £ evidently mea- 



Plugs in rough capillary tubes 11 
sures lengths in units comparable to the remobilization length Loq. We denote the plug 

length L p in these units as £ p . So in this section we will study long plugs with £ p ^> 1. 



As we approach the thin film region, 77 becomes arbitrarily close to 1: \r) — 1| -C 1. So 
we can simplify Eqn. (|2.5p as follows: d 3 r]/d^ 3 ~ 77 — 1. This autonomous equation has 
solutions of the form: r\ — 1 = cxp(i/£), where ^ 3 = 1. It has three roots v\ — 1 and 
z/2.3 = —1/2 ± \/3i/2, which correspond to three solution modes given below: 

^i(C)=exp(0, 

V>a(0=exp(-£/2) cos (v^/2) , 

Mt) = exp (-e/2) sin (V^/2) . (2.6) 

So general solutions are linear combinations of these modes, and we have: 

3 

, = l + ^^,(£), (2.7) 

where coefficients are determined from boundary conditions away from the thin film 
region. We note for future reference that for large positive £ , 77 is dominated by the 
growing ipx function. 



Far from the thin film region treated above, the following relation is valid: 77 = h/h* 3> 
1. For ve ry small Ca, lu brication approximation (\dh/dz -C 1|) is still applicable in this 



regime I (|Teletzke 19831 ). So Eqn. (|2.5|) can be approximated as follows: d 3 ij/d^ 3 ~ 0. 



Integration of this equation yields parabolic profiles: 



' 2 s 



lF,R 



F.R 



(2.8) 



where sup erscripts F and R stand for front and rear solutions, respectively. As defined by 
Teletzke| Jxeletzke 1983 ). axial (z-direction) positions of A and B in Fig(5J correspond 
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to axial positions of the lowest points of front and rear parabolas. And he defined the 
distance between A and B as the plug length L p . We will keep Teletzke's definition in 
the following discussion. 



Since both ends of the plug are spherical caps for small Ca, mean curvatures far 
from the thin film region must match those of spherical caps. As a result, the following 
constraint needs to be satisfied near the spherical cap: d 2 h/dz 2 = 1. We may rewrite this 
constraint in terms of reduced variables: (3Ca) 2 / 3 (h*)^ 1 ■ d 2 r]/d£ i ' 2 \^^ oo = 1. By using 
asymptotic solutions of rj given in Eqn. (|2.8p . we get: 

(^!.p^ = 1 . (2 . 9) 

h* y ' 

Since h* is fixed for a given plug, the above equation implies that P F = P R . So we can 
drop the superscript from now on. For a semi-infinite plug, we shall use Bretherton's 
result for h^. This leads to the limit value of P as the plug length goes to infinity: 
Pea = 0.643. As derived in detail in the Appendix of this paper, the dependence of 
the curvature P on the plug length £ p for finite long plugs is given by: P(£ p ) = -Poo — 
1.38 ^(^p) + 0.48(53(£ p ), where functions 82 and ^3 are linear combinations of modes ■02 
and ip3, defined in Eqn. (|A8p . By using the constant curvature constraint derived above, 
we have: 

h* = Pfo) 

hoo Poo 

= 1- 2.15 + 0.75 (2.10) 

This is the dependence of h* on the dimensionless plug length £ p . It can be easily trans- 
formed to represent the dependence on L p via the scaling relation: L p = £ p P(S,p) ■ 
(SCa) 1 / 3 . Thus the thickness h* oscillates around with an amplitude that decays as 
L p increases, and with a decay length of order L. Strictly speaking, the above derivation 
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Figure 4. Dependence of h* on the plug length L v . Solid line represents our data obtained via 
a perturbation analysis. Points are data obtained via Teletzke's method IjTeletzke 1983T ). 



is only valid in the regime: £ p 3> 1, where 82 and £3 are arbitrarily small. For illustrative 
purpose, we extrapolate our analytic results to Teletzke's short-plug regime and make a 
comparison between data obtained via two methods in Fig|4j Deviation between two sets 
of data is obvious for short plugs, where Teletzke's method is more accurate. For short 
plugs in the Teletzke's regime, contributions from all modes are important. Different 
modes interact nonlinear ly, which leads to a much stronger finite- length effect, as shown 
in Fig|3] Around L p ~ 6Loo the data obtained via Teletzke's method begin to merge 
with our curve. This is the regime where the prediction in Eqn. (|2.10p becomes valid. 



In this section, we completed the previously unstudied part of the finite plug speed 
problem. Together with Teletzke's method, we can now predict how fast plugs of different 
length will move in smooth cylindrical tubes. Evidently the speed dependence on the plug 
length is governed by the solution modes V'1,2,3 in the thin film region. For long plugs 
studied in this section, their speeds are the most sensitive to the mode with the longest 
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Figure 5. A). A sketch of a capillary tube with azimuthal roughness. B). A plug in the rough 
tube showing the quantities e(0), h(z,9), qe, L* defined in the text. 

wavelength: X max = max(|i?e(^ 1 )|). In smooth tubes, this X max is much smaller than 
the tube radius R for small capillary numbers of interest. The asymptotic corrections of 
the finite plug speed is exponentially weak, with a decay length X max . For very long plugs 
with L p >• Xmax, the speed is almost the same as that of the semi- infinite one. As dis- 
cussed in the next section, tube roughness can modify the spectrum of solution modes in 
the thin film region dramatically so that X max may become comparable to or larger than 
R. This is crucial in understanding some interesting experimental phenomena involving 
the relative motion between finite plugs of different length in non-smooth capillary tubes. 



3. Relaxation modes in the thin film region with tube roughness 

In this section, we investigate the shape of the thin film region of semi-infinite plugs in 
tubes with roughness. Here tube roughness is defined as surface deformation of the tube 
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wall, whose shape is not cylindrical. For simplicity, the only type of roughness that we will 

consider in the following discussion is azimuthal roughness. A tube with the azimuthal 

roughness preserves its translation symmetry in the axial direction, but cross sections are 

no longer circular. A sketch of such a tube is plotted in FigJSJA). The cross section of the 

rough tube is parameterized in the polar coordinate system as: R{8) = R + e(8), where R 

is the smooth tube radius and 8 is the polar angle. The roughness function e(8) measures 

the deformation. In order to isolate the regime of small roughness we may express e(8) as: 

e{8) = ea(8). Here e is some small parameter much less than the Bretherton's thin film 

thickness h^, and a{8) is a function of order unity. The tube roughness breaks rotational 

symmetry in the azimuthal plane. As a result, the plug shape now depends on both the z 

and the 8 coordinate: h(z,8). Local lubricating film thickness is defined by the distance 

between the tube wall and the plug surface: h(z,8) + e(8), as shown in Fig[5]jB). As in 

the previous section, we nondimensionalize the system by using the smooth tube radius 

R, the plug speed U and the capillary pressure j/R- As derived later in this section, the 

semi-infinite plug shape in the thin film region can be expressed as: 

h(z,8)~H 00 +J2CkM6)expM (3.1) 

k 

where ipk and Vk are eigenfunctions and eigenvalues governing the relaxation of the plug 
shape to the cylindrical asymptote . Relaxation mode coefficients C'k are determined 
by boundary conditions away from the thin film region. In Section 4, we will estimate the 
order of magnitude of Ck near the transition between the thin film and the remobilization 
region. 
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3.1. The plug shape equation in rough tubes 



As in the case of the smooth tube, the plug shape is dictated by the local carrier fluid 
flux q\z, 9): 



q(z,6) 



h(z,8) 



v(z, 6, r) dr , 



(3.2) 



-e(8) 



where v(z, 6, r) is the carrier fluid velocity. As shown in FigJSjB), if has both an axial 
component q z and an azimuthal component qg. In the limit of lubrication approximation, 
a nd with the same boundary conditio ns as in the smooth case, these fluxes are given by 



(jDeGennes. Brochard fc Quere 20031) 



Qz 
<lo 



(h + e)-^±^Q z P, 



3Ca 



(h + ef 
SCa 



d e P, 



(3.3) 



where P is the capillary pressure on the plug surface. In our reduced units (jHowelll (|2003h ): 
P = —h — (dl + dg)h. Eqn.(|3.3p resembles t he usual Poiseuille's law for thin layer fluid 



fluxes | (jDeGennes. Brochard fc Quere 20031 ). Rearranging terms in the equation of q z 
we get: 



3Ca 



(q 2 + h + e)~d z {l + d 2 e )h. 



(3.4) 



(h + e)3 

This is the counterpart of the Teletzke's plug shape equation, Eqn. (|2.4| ). in rough tubes. 
Without roughness, q z is a constant and this equation is sufficient to determine the plug 
shape h(z). However, h now depends on both 6 and z coordinates, the flux q z is no longer 
constant and another equation describing the variation of q z is needed to determine the 
plug shape. This additional equation is the incompressibility constraint: d z q z +dgqg = 0, 
which leads to the differential equation: 



o z q z = de —ttz — doP 



3Ca 
(h + ef 
3Ca 



d e [(1 + dj) h + 2 z h] 



(3.5) 
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Combining Eqn. (|3.4[ ) and (|3.5|) . we get a complete set of differential equations dictating 



both the plug shape h and the flux q z 



d z h 


= b! 


d z ti 


= h" 


d z h" 


= D- 


d z q z 


= -dt 



(3.6) 



The first two equations are simply definition equations of variables h! and h" . And the 
parameter D is given by: D = (h + e) 3 /(3Ca). Evidently, z and 9 derivatives are sepa- 
rated to two sides of the equation. 



In the smooth tube the plug surface has rotational symmetry, so the last equation of 
Eqn. (|3.6|) is simplified as: d z q z — 0. Obviously, it implies a constant flux q z . Defining this 
flux as — h* and using it in the third equation of Eqn. (|3.6[) . we get: 

„ 3 3Ca ■ (h~h*) 

d z h = x — - d z h . (3.7) 

Finally, we saw in the previous section that in the region of interest h varies over a length 
scale of about L ~ Ca 1 / 3 . So the slope term d z h in Ean. pjp is negligible compared to 
the other two terms for small capillary numbers. And the Teletzke's smooth plug shape 
equation, Eqn. (|2.4l) . is obtained after dropping d z h. 



Returning to the case of nonzero e(9), for convenience, we define a plug-shape vector: 
ip = {h, b! ', h", q z + e} T , where the superscript T indicates the matrix transpose oper- 
ation. Then we can rewrite the above set of differential equations in a concise matrix 
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d z ip = 



( 




D- 1 

y -d e Ddg{l + d 2 e ) 

= M(h)ip. 



1 


<l + d$) 






1 

D- 1 
-d e Dd e 



v 



(3.8) 



The nonlincarity of this equation arises because of the dependence of D on h. In the 
following discussion, the roughness function e(0) is assumed to be of a simple sinusoidal 
form: 

e{6) = e cos (m6) , (3.9) 

where m are positive integers indicating different roughness modes. The roughness am- 
plitude eo is assumed to be small in the regime of interest: eo hoo- As discussed in 
detail later, Eqn. (l3.8p can be reduced to a linear form in the thin film region. More 
general roughness functions can be expanded in terms of the Fourier series and analyzed 
correspondingly. 



3.2. Asymptotic shape of semi-infinite plugs 

In this section, we consider the asymptotic shape of the front region of a semi-infinite 
plug of arbitrarily long length L p ^> 1 . A similar analysis can easily be applied to the rear 

region of the semi- infinite plug. One immediate concern wit h these long plugs is that they 

may be unstable against fission via the Rayleigh instability! (jDeGennes. Brochard fc Quere 20031 ). 
This instability is independent of the steady state effect that we discuss here. Indeed, for 
small capillary numbers, the growth rate of the Rayleigh instability becomes negligible | 
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( DeGennes. Brochard fc Quere 20031 ). 
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As with the smooth tube case discussed above, we may fix the origin of z coordinate 
at the transition between the thin film and the remobilization region. Far inside the thin 
film region where z <C — 1 , the front semi- infinite plug restores the translation symmetry 
in the z direction. It implies vanishing z derivatives in this limit, i.e. d z ip(oo) — 0. Evi- 
dently, h'(oo) and h"(oo) are zero. The third equation of Eqn. (|3.8|) leads to the condition: 
h = —q z — e. And the asymptotic plug shape h(po) is determined by the last equation 
of Eqn. (|3.8p : —dgDdg (1 + dg)h = 0. One obvious solution is that h(oo) = H^, where 
iJoo is a constant. This solution corresponds to a cylindrical plug centered inside the 
tube. There are two other solutions: Hi cos 9 and H2 sin#, which correspond to cylindri- 
cal plugs shifted laterally inside the tube. The translation symmetry far inside the thin 
film region leads to the cylindrical plug shape with constant capillary pressure on the 
plug surface and vanishing flux qg. For convenience, we focus on the centered asymptote. 
The implication of non-centered solutions will be addressed in the Discussion section. 

The centered solution Hc^ measures the magnitude of average q z per unit circumfer- 
ence: 



So -f/oo has the same physical meaning as h x in the smooth tube case, and dictates the 
relative motion between the plug and the carrier fluid: U/V = 1 + 2Hoo. For tubes with 
small roughness, the difference between and hoc is also very small: \Hoo — hoo\ *C /loo- 




= -H t 



OO 



(3.10) 
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The leading order correction in Hoo due to the tube roughness appears in the second or- 
der of e. 



3.3. Eigenmodes of relaxation in the thin film region 

At finite z values, the semi-infinite plug shape deviates from its cylindrical asymptote. 
The plug-shape vector ip is rewritten in the following form: 



ip = ip(oo) + Sip 
= {Hoo + Sh, Sh', Sh", -Hoc + Sq z } q 



(3.11) 



where Sip = {Sh, Sh', Sh" , Sq z } measures the amount of shape correction. In the thin 
film region, Sip satisfies the relation: | Sip 1 | << | ip l {oo) | , where the superscript i indicates 
different vector components. So the plug shape equation (|3.8|) can be linearized by using 



the approximation: D ~ Di = (Hrx, + e) 3 /(3Ca), and we get: 



d z Sip 



1 

1 

D- 1 -(l + 9 e 2 ) 

^ n -d B D 1 d g 



^r 1 



y-dgD.deil + d^) 



Sip 



(3.12) 



Evidently, we can solve the above equation via the technique of separation of variables. 
The general solution of Sip is given by: 



Sip = ^2C k ip k {6) exp (v k z) , 



(3.13) 



where ipk{&) and z//. are eigenvectors and eigenvalues of the matrix operator Mi, with 
M-iipk(9) = Vkipkiff). Thus in the thin film region, ipk(6) exp(^fcz) are the eigenmodes 
dictating the relaxation of the plug surface to its cylindrical asymptote. The speed of 
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Figure 6. Zeroth harmonic spectrum of the semi- infinite plug shape in the thin film region for 
tubes with small sinusoidal roughness. Capillary number's value: Ca — 10 -3 . 

relaxation is determined by the decay length of each mode Afc, defined as Afc = \Re(l/i/k)\- 
The larger the value of Afc, the slower of the relaxation process. Relaxation modes with 
positive i?e(l/Vfc) are dominant for the front semi-infinite plug. Likewise, negative ones 
are dominant for the rear semi-infinite plug. 



In the regime of vanishingly small roughness e((9), Eqn. (|3.12p can be simplified further 
as follows: 

/ 

10 

10 



d z Stp 



Do 1 -a + d 2 e ) Do 1 

y -D (d 9 S + d%) -D d 2 e 



5 th 



(3.14) 



where the coefficient Dq is now given by: -Do = ft^/ (3Ca). Evidently the matrix operator 
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Mo has no explicit dependence on the 6 coordinate, so its eigenvectors correspond to dif- 
ferent harmonics in the Fourier series of 9. Due to the form of the roughness function e{6) 
chosen, we may focus on cosine series solutions with the same symmetry: ^ C n cos (nmff). 
Here n is the harmonic number. Up to the first order in the roughness e, we need to an- 
alyze the zeroth and first harmonics, with each harmonic leading to four eigenmodes: 



5$ ~ (5?A (0) +5%f {1) 

= E C k 0) ^ exp (^ 0) z) + E cos (m6) exp (u^z) , (3.15) 



fe=0 



k=l 



where superscripts (0) and (1) stand for the zeroth and the first harmonics respectively. 
The zeroth harmonic solution 5ip^ preserves the rotational symmetry and satisfies the 
equation: 

/ 

10 
10 



\ 



-1 D. 



(I 



'0 



The above matrix operator is a four-by-four matrix and does not depend on the form of 



V 



6il> 



(0) 



(3.16) 



/ 



the roughness function. So the zeroth harmonic spectra of different roughness functions 
are the same. For small capillary numbers, four eigenvalues of this matrix are given 
by: „C°> = 0, = D- 1/3 , 4 0) = J D - 1/3 (-l/2 + * v / 3/2), and = ^ 1/3 (-l/2 - 
zn/3/2). There is no plug surface relaxation associated with the mode v ^ . It corresponds 
to the degree of freedom to fix the radius of the cylindrical asymptote as z — > oo. The 
other three modes , and correspond to the ones found in the previous section 
for the smooth tube, up to a difference in the choice of units. For illustrative purpose, the 



zeroth harmonic spectrum is plotted in Fig|51 Obviously, and v£' > are the relaxation 
modes with the longest decay length AmL in this spectrum. 



,(0) 
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Figure 7. First harmonic spectrum of the semi-infinite plug shape in the thin film region for 
tubes with small sinusoidal roughness. Cases with m =5, 6, 7 are plotted. Capillary number's 
value: Ca = 10" 3 . 

A major difference between the rough tube and the smooth one comes from the first 
harmonic solution Sif)^. It satisfies the equation: 



d z 6$ w = 



Do 1 



Dq (TO 2 — TO 4 ) 



1 



-1 + m 2 





1 


Dq TO 2 






Do 




5^ 



Four eigenvalues are solutions of the characteristic equation: 

i/ 4 + (1 - 2to 2 ) v 2 - D^ 1 v + to 4 - to 2 = . 



(3.17) 



(3.18) 



The first harmonic spectra with the roughness mode number to =5, 6, and 7 are plotted 
in FiglT] Compared to the zeroth harmonic spectrum shown in Figj6l the mode with the 
longest decay length AmL now lies near the origin and has much longer decay length. For 



small values of m, i.e. long wavelength roughness, AmL becomes comparable to or even 
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larger than the tube radius R. Obviously, Xmax > ^max in the regime of interest. So XmLx 
is the new slowest decay mode in rough tubes after analyzing both the zeroth and the hrst 
harmonic spectra. So we will drop the superscript (1) and refer to it as X max from now on. 



To better understand the origin of these slow relaxation modes in tubes with rough- 
ness, we propose a following physical explanation of X max . The remobilization region 
ends at about distance z ~ from the spherical cap. At this transition between the 
remobilization region and the thin film region, the plug shape has not reached its cylin- 
drical asymptote yet. We assume that the plug surface has a typical height variation of 
about Ah. A sketch of such a plug cross section is given in FigJSjB). Non-cylindrical plug 
shape leads to unbalanced capillary pressure P in the azimuthal plane and a non-zero 
lateral flow qg. In the thin film region, the plug surface relaxes toward the cylindrical 
asymptote via this flow qg. Suppose that the relaxation length scale is given by L* . It can 
be estimated as the product between the carrier fluid speed U in this region and a relax- 
ation time scale T over which qg smoothes the plug surface. In the following analysis, we 
will estimate the length scale L* and obtain its scalings with other system parameters. 
For this part of discussion only, we will work in the lab units. 



The expression of the lateral flux qg is given by Eqn. (|3.3p . In the lab units, we have: 

As mentioned above, the capillary pressure P is no longer a constant in the azimuthal 
plane for the non-cylindrical plug shape. The length scale in the 9 direction is de- 
fined by the roughness wavelength lg ~ R/m. So we may estimate 9 derivatives as: 
i? -1 dgh ~ (m/R) - Ah. Since the z-direction plug surface adjustment mostly happens in 
the remobilization region, we expect that 9 derivatives are typically much larger than z 
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derivatives in the thin him region. The limit of this assumption will be derived later. So 



the capillary pressure P is approximated as follows 

P = - 7 



R 2 

7 7^- (3-20) 



Ah 

(R/m) 



By dehnition, the surface relaxation time scale T is given by: 

m _ Ah M (i?/m) 4 



(3.21) 



R 1 d e qe ih^ 

The relaxation length scale L* is estimated as the product between U and T: L* ~ 
Ca R 4 / (h^ m 4 ). We may simplify this estimation of L* further by using the scaling of 
hoo and get: 

11 (3.22) 



Ca ■ m 4 

For a self-consistent derivation, m's value is constrained by the approximation used in 
Eqn. (|3.20p . This is equivalent to the condition: R^ 1 dgh 3> d z h, which yields the relation: 
lg CL*. By using the estimation of L* obtained above, we get the valid regime of the 
roughness mode number: m -C Ca -1 / 3 . It corresponds to capillary tubes with long wave- 
length roughness. Obviously L* is much longer than the smooth remobilization length 
Lqo in this regime of roughness. 



The length scale L* is actually the longest decay length X max in the eigenvalue spec- 
trum derived previously. It corresponds to the smallest solution of Eqn. (13.181) , which goes 
to zero as Ca approaches zero. For this root, the D^ 1 v term dominates the higher power 
terms of v, and Eqn. p.l8[) simplifies as follows: 

-Do 1 'X^+m i -m 2 = 0. (3.23) 

For an order of magnitude estimation, we may drop the term m 2 relative to to 4 and get: 
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A mn i ~ Ca/lh^ m 4 ). This is the same scaling relation as that of L* , up to a difference 
in the choice of units. Thus L* and A max measure the same length scale in rough tubes. 
The scaling argument presented above physically explains the origin of large A mQX values 
in the rough tube. This new length scale reflects the slow relaxation of the plug shape via 
a vanishingly small lateral fluid flux qg. Historic ally, people ha s observed some evidence 



of this new length scale in non-cylindrical tubes 



Wong 1 



l|l995h . 



4. Amplitude of the plug surface bumpiness in the remobilization 
region 

In the previous section, we analyzed the semi-infinite plug shape in the thin film region 
and found that the plug surface decays to its cylindrical asymptote via different relaxation 
modes. To proceed further, we must understand how the rough tube wall imposes a non- 
cylindrical plug shape in the first place. For small roughness, the plug shape is dominated 
by the zeroth and the first harmonic solutions: 

3 4 

h(z, 6)^H OG +Y, C£ 0) exp (^ 0) z) + C£ ] cos (m0) exp z) , (4.1) 

fe=l k=l 

where the coefficients C^ ' and are determined from the boundary conditions away 
from the thin film region. In this section, we will treat the front and the rear semi-infinite 
plugs simultaneously, with the understanding that only the modes with the correct sign 
of Re(l/uk) are excited for a single semi-infinite plug. The z coordinate origin is fixed 
at the transition between the thin film and the remobilization region as in Section 2. In 
that section we have shown that the zeroth harmonic coefficients are of the order of 
/loo. We will see later in this section that the first harmonic coefficients C^p are of the 
order of eo. The exact numerical solutions will be presented in our future work. 
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In the remobilization region, the governing equation of the plug shape h(z, 9) is derived 

by combining the last two equations in Eqn. (|3.6p : 

d z (h + e) - A- [d z (h + efd z + d e (h + efdg] (1 + 8 2 z + d 2 )h = . (4.2) 

This equation can be simplified further by noticing the fact that in the remobilization 
region, h varies over about (SCa) 1 / 3 in the z direction and about mT 1 in the 9 direction. 
For long wavelength roughness with m <C (3Ca) -1 / 3 , the following estimation of deriva- 
tives is valid : d z h dgh. In other words, in the remobilization region the plug shape 
h(z, 9) varies much more rapidly in the z direction than in the 9 direction. So we may 
keep only z derivatives in Eqn. (|4.2l) . and get: 

d z h=^-[d z {h + efdl]h. (4.3) 

Bretherton's smooth tube solution ho(z) satisfies the above differential equation with 
e = 0. In the regime of small roughness, Eqn. (|4.3p can be solved perturbatively by expand- 
ing around the smooth tube solution ho(z): h(z,9) = ho(z) + u^°\z) + u^ x \z) cos(m9). 
Here u' ' and are the leading order plug shape correction attributed to the tube 
roughness. For vanishingly small roughness, the relation \u( 0,1 ^ /h$\ <§; 1 holds for arbi- 
trary z. Since we are interested in the coefficients of the first harmonic solution, we will 
focus on the shape correction uW for the rest of this section. Keeping only the lowest 
order approximation, satisfies the following equation: 

(d z hl dl) u (1) + (3h 2 d 3 z h - 3Ca) d z u (1) + [d z (3h 2 d 3 z h )] {u (1) + e ) = . (4.4) 

This is an inhomogeneous differential equation. The correct boundary conditions of 
are that it goes to zero as z approaches the spherical cap, and matches the relaxation 
modes in Eqn. (|4.ip as z approaches the thin film region. The general solution of Eqn. (|4.4l) 
takes the form: u^(z) — eo (f(z) — 1). Here f(z) is a dimensionless function satisfying 
the homogeneous equation and goes to 1 as z approaches the spherical cap. By using the 
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rescalings as in the smooth tube: ho = hoo ■ 770, and z — ■ £, we get the following 
equation of /: 



After the rescaling both 770 and its £ derivatives are of the order of unity in the remobi- 
lization region, so all coefficients in Eqn. (|4.5l) are of the order of unity. Together with the 
boundary condition of / near the spherical cap, the solution f(z) must be of the order 
of unity in the remobilization region. By matching with relaxation modes in the thin 
film region, we get an order of magnitude estimation of the first harmonic coefficients: 
ci ~ eo- Thus the semi-infinite plug shape in the thin film region can be rewritten as 
follows: 



where all coefficients C k and are of the order of unity. The above estimation is 
valid for both the front and the rear semi-infinite plug. 

5. Finite plug motion in rough tubes 

In this section, we qualitatively study the dependence of the finite plug motion on 
the plug length in the rough tube, following the logic of Section 2. For the finite plugs, 
relaxation modes of both signs of Re(l/uk) are excited in the thin film region. Far inside 
the thin film region where the plug shape is almost cylindrical, all the other relaxation 
modes are vanishingly small compared to the slowest decay mode X m ax • Thus we may 
focus on the effect of \ ma x- By defining a z coordinate whose origin is at some position 
far inside the thin film region, the plug shape is approximated as follows: 



a|) / + (3 % 2 df % - 3) dtf + ft(3»7§ df %)] / = . 



(4.5) 




(4.6) 




(5.1) 
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where z = z + £. The coefficient C is of the order of unity as derived in the previous 

section. Here the variable £ measures the half length of the finite plug. For plugs used in 

microfluidic experi ments, £ is usually comparable to or larger than the tube radius R \ 



(jlsmagilov fc Yina ). As mentioned previously, the relative motion between the plug and 
the carrier fluid is governed by the azimuthal average of q z : 

2tt 



u/v=i~- r 

I" JO 



q z d9 



= l-2q. (5.2) 

Far inside the thin film region, the plug shape is almost cylindrical, which leads to the fol- 
lowing estimation of the carrier fluid fluxes: q z ~ — h—e, and qg ~ 0. By using the approxi- 
mate plug shape in Eqn.(dTT]), we have: q z (z, 9) ~ —H 00 — Ce{9) exp[(z — £)/ \ max ] — e(9). 
Evidently, the average flux q vanishes in this linear approximation. Thus any effect on q 
appears in the second or higher order of the roughness: q + Hoo — B e\ exp(— 2l/A max ) + 
O(efj). The asymptotic dependence of the plug speed on the plug length is exponentially 
small for long plugs, and we have: 



U/V ~ 1 + 2H a 



Be 2 

1 + -tf^ exp(-2£/A max ) 

-Hoo 



(5.3) 



The above form of length dependence implies an interesting resonance effect in the rough 
tube. To illustrate this effect, we may estimate the speed difference between two plugs 
with length l\ and 1-2,. We also assume that the length difference between these two plugs, 
defined as A£ = £2 — £1, is small such that |A£/€i| <§; 1. By using Eqn. (|5.3p . we derive 
their speed difference as follows: 



(U 2 - U{)/V ~ 2Be 2 [exp(-24/A ma:c ) - exp(-2l 1 /A moa )] 

~ 2Be 2 (-^-) ■ [exp(-2^/A max ) • (2£ 1 /X max )} , (5.4) 
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The function of A max in the last square bracket above is not monotonic. It has a max- 
imum at X m ax = For plugs used in microfluidics measurements, we have: l\ ~ R. 
Thus the speed difference of these plugs is maximized in rough tubes with X ma x ~ 2R. By 
using the scaling of X max derived in Section 3, we get the roughness mode of resonance: 

In the rest of this section, we will use an example to illustrate how the tube roughness 
enhances the sensitivity of the plug motion to the plug length. Suppose that we have 
two finite plugs of length: l\ = 0.5R and £ 2 = R- Strictly speaking, Eqn. (|5.3p is valid for 
the regime of vanishingly small roughness. However, we expect this estimation to be at 
least qualitatively correct for tubes with large roughness Thus for computation simplic- 
ity, we may set the coefficient B eg/i/oo to be 1. In the smooth tube case, ffoo = hoo and 
Xmax = 2Loo. And we get the speed difference between these two plugs less than 0.01% 
with Ca — 0.001. Thus the finite plug length effect is negligible in the smooth tube. On 
the other hand for tubes with long wavelength roughness, as shown in Fig|6l we have 
parameter values: ~ and \ ma x — R- The plug speed difference can be up to a 
few percent, which is significantly larger than that in the smooth tube. 

6. Discussion 

In the previous sections we have argued that roughness in a microchannel can quali- 
tatively alter capillary motion in the channel. We demonstrated the new effects in the 
context of an air plug being pushed through the channel. Roughness creates new modes 
of relaxation of the fluid interface. These modes make major changes in the Bretherton- 
Teletzke mechanism that sets the film thickness and governs the plug's speed. In partic- 
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ular, they alter the nature of the the remobilization region that dictates the thickness of 

the liquid film between the plug and the wall. These modes can thus alter the effect of 

other variations in the system, such as the addition of surfactants or the replacement of 

the air plug by a viscous fluid. 



Our focus above has been to consider the effect of plug length on its speed. This basic 
effect has practical consequences as noted in the Introduction. Microfluidic devices for 
chemical processing create long trains of droplets that need to remain separated over long 
distances. However in practice these droplets can have slight differences of speed. This 
can make droplets merge, thus disrupting the desired behavior. Since the length of these 
droplets is not precisely controlled, the length differences could in principle be respon- 
sible for the undesired speed differences. However, such a length effect is incompatible 
with the classic Bretherton explanation for the speed. We showed above that roughness 
changes this conclusion qualitatively. 



This effect complements previously-considere d effects of channel shape. It d iffers from 
the roughness effect identified by Homsy et al. | ( Krechetnikov fc Homsv 2005 ). in which 
roughness alters the effective boundary condition at the wall. We showed here that rough- 
ness can also have long-range effects on the thin film behavior. These effect s do not requir e 



the qua l itativ e changes in channel shape found, e.g., for square channels (jWong 1 



Wong 2 



dl995h : 



(|1995T )V Remarkably, even high-wavenumber roughness can alter the long-range 



effects if the capillary number is sufficiently small, as shown in Section 5. 



Roughness may be viewed as a relevant perturbation to the plug viewed as a dynamical 
system. The fast relaxation to the thin film within a narrow remobilization region seen 
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in a Bretherton plug depends crucially on its assumed azimuthal symmetry. When this 
symmetry is broken by roughness, the relaxation takes on a qualitatively new behavior. 
In this sense the symmetric plug is unstable against the generic un-symmetric behavior 
discussed above. For this reason it appears important to consider roughness when esti- 
mating any new aspect of plug behavior. 

Though our work suggests that roughness is important, our quantitative exploration 
of this importance has been very limited. First, we considered only azimuthal roughness, 
ignoring any z dependence of the wall position. We believe the nonlocal effects identified 
above are not qualitatively altered for generic roughness, we have not addressed this issue 
explicitly. We consider only perturbativc effects of roughness in lowest order. Thus we 
have no systematic calculation for the speed, which requires the next order of approx- 
imation. Our analysis of the remobilization region was limited to a scaling discussion, 
with no systematic results. We didn't carry out more detailed calculations because the 
perturbation approach itself has serious flaws. It becomes qualitatively inadequate for 
low capillary number, where the unperturbed film thickness becomes smaller than any 
roughness. In this regime the interface is mainly supported by asperities where the film 
thickness is much smaller than its average. In order to get quantitative information on 
how roughness and length affect speed, we plan a numerical integration of Eqn. (|3.8[) . 

Our analysis unearthed a peculiar behavior for the lowest-wavenumber perturbation, 
namely m = 1. This mode amounts to a lateral displacement of the plug in a circu- 
lar channel. Remarkably this perturbation does not relax since it does not alter the 
circular shape of the plug. Centering effects of particles in tubes have been discussed | 
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( Brenner 19631 ). but we are not aware of an analysis suited to Bretherton plugs. 



The roughness effects shown here suggest a new regime of interest in microfluidic capil- 
lary flow. In this rough-channel regime, the symmetry of the Bretherton plugs is lost, but 
the fluid coating the interface is still a thin film in lubrication flow. Thus th ese channels 



(19951)), where the 



are distinct from, e.g., the square channels discussed by Wong (jWong 1 
lubrication approximation is inadequate. In addition to capillary number and the plug 
aspect ratio, two dimensionless variables characterize a system in this regime. The first is 
the roughness wavelength R/m relative to the Bretherton film thickness h^. The second 
is the amplitude eo relative to this thickness. We have noted above the special effects to 
be explored when the decay length is comparable to the plug length. We have also sug- 
gested an asperity-dominated regime to be expected when the dimensionless amplitude 
becomes large. Naturally the chief experimental impact of these effects comes when the 
plug is replaced by a droplet of viscous liquid. The effects anticipated above should all 
have counterparts for such liquid droplets. 



Our predictions, though qualitative, are subject to experimental tests. An obvious test 
is to create tubes with controlled corrugated profiles and verify that when the roughness 
amplitude is comparable to the predicted Bretherton film thickness, the speed increment 
U—V changes significantly relative to that of a smooth tube with the same cross section. 
One can also create corrugations with a wavelength where sensitivity to plug length is 
expected and look for significant variations in U — V when the roughness amplitude is 
comparable to the Bretherton thickness. Such tests could be done at the microfluidic 
scale or at larger scales. 



([Ismagilov fc Yingf ). For such flows, the 
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More primitive observations using microfiuidic channels are also suggested by our anal- 
ysis. One may readily prepare circular tubes wit h 200 micron cross section and control- 
lable flows with capillary numbers of order f 0~ 5 
predicted Bretherton thickness is 60 nanometers, so that roughness amplitudes of this 
order should be significant. AFM measurements on these tubes show roughness with m 
of order 10 2 and amplitude of order 50 nanometers. This amplitude is evidently large on 
the scale of interest. Since the roughness varies randomly along the tube, one expects 
statistical fluctuations in the speed U along the tube. Some of these fluctuations can be 
due to differences in cross sectional area, so that they cause the fluid speed V to change 
with position. In addition one expects variations in U — V owing to our mechanism. The 
indication that these variations are caused by roughness is that U — V is determined by 
position. Thus e.g. large U — V values occur at particular places along the tube. 



7. Conclusion 

The subtle interplay of nonlocal forces that controls a microfiuidic droplet is a classic 
example of the power of simple capillary flow to produce nontrivial structures. It appears 
from our analysis that weak corrugation of the guiding channel can add unexpected rich- 
ness to this picture. The addition of transient effects, binary interaction between droplets 
and surfactants seems likely to add even more to this richness. Exploring these effects 
appears promising both for microfiuidic applications and for basic understanding of what 
capillary flow can do. 
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Appendix A. 

In this appendix, we derive the dependence of the rescaled curvature P on the finite 
plug length £ p . All quantities are in the dimensionless units defined by the rescalings in 
Section 2. To set the stage for our treatment of finite plugs, we analyze the behavior of 
rj for the semi-infinite plug treated by Bretherton. We proceed from a semi-infinite plug 
with the front cap. Far inside the thin- film region, |r/— 1| <C 1 and we may use the general 
solution in Eqn. (|2.7p . In the front region of the thin film, the contributions from ip2 and 
ip3 are arbitrarily small relative to that from tpi. So the plug shape takes the form: 

V = l+CiM0- (Al) 

Consistency with this form implies: rj^ — t](, = V ~ 1> where 77^ and r)£ stand for the 
second and the first derivative of 77 with respect to £ respectively. We may thus find 
the profile rj (£) by integrating forward from an arbitrary origin to £ — > +00, starting 
from an initial value of rj very close to 1. The asymptotic curvature rj^(+oo) is necessar- 
ily finite and independent of axial translation in the small-amplitude limit. Specifically, 
7/^ (+00) Poo = 0.643. It is convenient to define £p as the £ for which the ipi contri- 
bution in Eqn. IjA 1|) extrapolates to 1, i.e. Cii/ji(£,f) = 1- Thus the thin film region of 
Eqn. (|A II) is where £f — £ S> 1, and ?7(£) = 1 + e^^^ F . Then from the ?7(£) profile thus 
determined numerically, one finds that Teletzke's point B has position: £b — £f — 0.17. 
Evidently, rj^ (+00) is independent of £77. 

A similar analysis can be applied to a semi-infinite plug with the rear cap. Here, the 
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ipi contribution is negligible in the rear region of the thin film. So the plug shape takes 
the form: 

v = i + c 2 M0 + c 3 M0- (A 2) 

As with the front region, we may find the rear r\ profile by integrating backward from 
the thin-film region where r\ — > 1. Far inside the thin- film region the ip2 and ip3 con- 
tributions are necessarily small, but their ratio varies sinusoidally with position, with a 
wavelength comparable to the remobilization length L^. As above, we define a £r such 
that C^V^CCr) = 1- Then in the thin film region n takes the form: 

V = 1 + e-«-««>/ 2 • [cos (V3(£ - £ fl )/2) + C 3 sin (V3(£ - £ R )/2)] • (A 3) 

To determine 7y(£) in the nonlinear regime, we integrate towards negative £ starting from 
initial conditions compatible with this form. An initial value £0 is chosen such that 77 (£0) 
is very close to 1. Compatibility with Eqn. (|A3| ) dictates the values of ?7j(£o) and i]^(£,o)- 
Then we integrate towards £ — > —00 to obtain the profile. As with the front profile, the 
asymptotic rear curvature r/^ (—00) is necessarily finite and independent of translation 
in £. However, this r)^(-oo) depends on the choice of C3 values. The correct profile is 
that which matches the front curvature rj^ (—00) = r/^ (+00) = 0.643. This requirement 
fixes the value of C3 in Eqn. IjA 3j) . We find numerically that C3 = —0.85. As with the 
front profile, we may determine the rear Teletzke point £a by numerically solving the 
quadratic coefficients of Eqn. (|2.8[) and get: £,4 = £# + 3.37. 



We now focus on long plugs of finite length: £ p 3> 1 . To the lowest order approximation, 
the above determination of Teletzke points allows us to relate the dimensionless plug 
length £ p to £f — £r as follows: 

£f-6? = £ P + 0.17 + 3.37 
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= Z P + 3.54 . 
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(A 4) 



Obviously £p —£r^S> 1 in the regime of interest. For finite plugs, the thin film region con- 
tains nonzero contributions from all three modes ipi, ip2, and ip^. Using our conventions 
for the f^'s above, rj takes the following form in the front region of the thin film: 

rj = l + e^ F +e - {£ --t R)/2 ■ |cos(V3(£-£ R )/2) - 0.85 sin^^ - &)/2) 

= i + Mi-ip) + Mi-iR)- °- 85 Mi ~ &) 

= i + Mi - ip) + Mi -ip + ^F- in) - o.85 Mi -ip + dF- in) , (as) 

where higher order perturbations due to finite plug length are neglected. The admixture of 
ip2 and ip3 contributions alter the behavior at large £ and thus perturbs the asymptotic 
curvature P. In order to find this perturbation, we must first express the ^>'s using a 
common origin £p. By using the addition properties of sine and cosine functions, we 
have: 

Mi -iF + ir- in) = Mi - ip)Mip - in) - Mi - ip)Mip - in) 
Mi ~iF + i F - in) = Mi- ip)Mip -i R ) + Mi- i F )Mi F - in) . (a 6) 

Plugging the above expressions into Eqn. (|A5|) and using Eqn. (|A4|) to substitute £p — in, 
we get: 

ry = i + Mi - ip) + s 2 (Q ■ Mi - ip) + Hi P ) ■ Mi - ip) , 

where 82 and ^3 are functions of the plug length £ p defined as follows: 

V> 2 (£ P + 3.54) 
^ 3 (£ P + 3.54) 

In the regime of interest £ p ^ 1, 82 and 83 are very small: |£ 2 ,3(£p)| ^ 1- Infinite plugs 
correspond to the limit case where #2,3 (+00) — > 0. Similar to the above analysis, we 
integrate towards £ — > +00 from some point £0 far inside the thin film region with 
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initial conditions compatible with Eqn. (|A 7|) . The asymptotic front curvature 77^ (+00) 
is necessarily finite and independent of small-amplitude translation in £. However, this 
ry^^(+oo) depends on the plug length £ p . Moreover, it is evident from the functional form 
of rj in Eqn. (|A 7[) that n^(+oo) has dependence on £ p only through functions 62^) and 
^3(^p)- By defining P = n^(+oo), we have: 



~p(o,o) + <5 2 (£ p )- 



dP 



(0,0) 



dP 



(0,0) 



p. 



dP 



hit?) 



(0,0) 



dP 



(A 9) 



(0,0) 



where we have used the fact that 62 and ^3 take vanishingly small values in the regime of 
interest and kept only the lowest order approximation. Partial derivatives (dP/d6%)\ofi 
and (dP/d5s)\o,o ar e sensitivity factors that we can determine numerically as follows: 



4.38, 



(0,0) 



dP 

dP 
d6~ 3 



Plugging these sensitivity factors into Eqn. (|A9l) . we get: 



= 0.48 . 



(0,0) 



P(Zp) = Poo- 1.38 + 0.48 ^fe) . 



(A 10) 



(All) 



This is the dependence of the curvature P on the plug length £ p stated in the text. 



REFERENCES 
Bretherton, F. P. 1961, J. Fluid Mech., 10, 166. 
Ismagilov, R. & Liu, Y., Private Communication. 
Gunther, A. & Jensen, K. F. 2006, Lab Chip, 6, 1487. 
Stillwagon, L. E. & Larson, R. G. 1988, J. Appl. Phys. 63, 5251. 
Schwartz, L. W. & Weidner, D. E. 1995, J. Engng. Math., 29, 91. 



Plugs in rough capillary tubes 39 
Kalliadasis, S., Bielarz, C, & Homsy, G. M. 2000, Phys. Fluids, 12, 1889. 
Mazouchi, A. & Homsy, G. M. 2001, Phys. Fluids, 13, 2751. 
Howell, P. D. 2003, J. Engng. Math., 45, 283. 

Ajaev, V. S. & Homsy, G. M. 2006, Annu. Rev. Fluid. Mech., 38, 277. 

Chen, J. D. 1986, J. Coll. Inter]. Sci., 109(2), 341. 

Einzel, D., Panzer, P. & Liu, M. 1990, Phys. Rev. Lett, 64, 2269. 

Miksis, M. J. & Davis, S. H. 1994, J. Fluid Mech., 273, 125. 

de Gennes, P. G. 2002, Langmuir, 18, 3413. 

Krechetnikov, R. & Homsy, G. M. 2005, Phys. Fluids, 17, 102108. 
Wong, H., Radke, C. J. & Morris, S. 1995, J. Fluid Mech., 292, 71. 
Wong, H., Radke, C. J. & Morris, S. 1995, J. Fluid Mech., 292, 95. 
Hazel, A. L. & Heil, M. 2002, J. Fluid Mech., 470, 91. 
Wenzel, R. N. 1936, Ind. Eng. Chem., 28, 988. 
Herminghaus, S. 2000, Europhys. Lett, 52, 165. 

BlCO, J., Tordeux, C. & Quere, D. 2001, Europhys. Lett, 55, 214. 

Teletzke, G. F. 1983, Ph.D. Thesis, Univeristy of Minnesota. 

Hodges, S. R., Jensen, O. E. & Rallison, J. M. 2004, J. Fluid Mech., 501, 279. 

de Gennes, P. G., Brochard, F. & Quere, D. Capillarity and Wetting Phenomena 

(Springer- Verlag, New York, 2003). 
Brenner, H. 1963, Chem. Eng. Sci, 19, 519. 



